package scu.maqiang.possion;

import scu.maqiang.fes.BVPType;
import scu.maqiang.fes.FES3T41;
import scu.maqiang.mesh.Mesh3T4;
import scu.maqiang.mesh.ScalarFunc;
import scu.maqiang.numeric.Direct;
import scu.maqiang.numeric.NewIterSSolver;
import scu.maqiang.numeric.SRMatrix;

/**
 * 三维泊松方程解, 周期性边界条件
 */
public class PoissonPeriodic3DT4 {

	public static void main(String[] args) {
		Mesh3T4 mesh2D = new Mesh3T4().cube(80, 80, 80);
		mesh2D.scale(Math.PI * 2, Math.PI * 2, Math.PI * 2);
		ScalarFunc sourceFunc = (xy, label, param) -> Math.sin(xy[0] + 0.25 * Math.PI) * Math.cos(xy[1] + 0.25 * Math.PI)
				* Math.cos(xy[2] + 0.25 * Math.PI);
		FES3T41 fs = new FES3T41(mesh2D);
		SRMatrix A = new SRMatrix(fs.GetNdof());
		double[] RHS = new double[fs.GetNdof()];
		fs.assembleStiff(new double[] {1.0}, BVPType.COMMON, A);
		fs.assembleSource(new ScalarFunc[] {sourceFunc}, null, BVPType.COMMON, RHS);
		double[] x = new double[fs.GetNdof()];
		//A.sort();
		//System.out.println(A.display());
		System.out.println("Apply periodic condition...");
		fs.applyBC_Periodic(A, RHS, Direct.Z, 1, 2);
		System.out.println("Apply periodic condition...");
		fs.applyBC_Periodic(A, RHS, Direct.X, 3, 4);
		System.out.println("Apply periodic condition...");
		fs.applyBC_Periodic(A, RHS, Direct.Y, 5, 6);
		//fs.applyBC_Periodic(A, RHS, new int[][] {{2, 4}, {1, 3}}, new Direct[] {Direct.X, Direct.Y});
//		System.out.println(A.display());
		//fs.applyBC_MBN(A, RHS, Direct.All, 0.0, 1, 3);
		//A.sort();
		//System.out.println(A.display());
		//System.out.println(A.isSymmetric());
		NewIterSSolver solver = new NewIterSSolver(A);
		solver.PCGSSOR(RHS, x, 1.5, 2);
//		DirectSSolver solver = new DirectSSolver(A, fs.getnBand());
//		NewIterSSolver solver = new NewIterSSolver(A);
//		solver.PLUSolve(RHS, x, 1);
		mesh2D.toTecplot("PeriodicSolution3D.dat", x);
	}

}
